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Sedimentation of the neutron rich isotope 22 Ne may be an important source of gravitational 
energy during the cooling of white dwarf stars. This depends on the diffusion constant for 22 Ne in 
strongly coupled plasma mixtures. We calculate self-diffusion constants Di from molecular dynamics 
simulations of carbon, oxygen, and neon mixtures. We find that Di in a mixture does not differ 
greatly from earlier one component plasma results. For strong coupling (coulomb parameter F" > 
few), Di has a modest dependence on the charge Zi of the ion species, Di oc Z~ 2 ^ 3 . However Di 
depends more strongly on Zi for weak coupling (smaller T). We conclude that the self-diffusion 
constant D^e for 22 Ne in carbon, oxygen, and neon plasma mixtures is accurately known so that 
uncertainties in Dn 6 should be unimportant for simulations of white dwarf cooling. 
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I. INTRODUCTION 

Observations of cooling white dwarf (WD) stars pro- 
vide important ways to date stellar systems [I]. White 
dwarf cooling involves simple physics because the energy 
from nuclear reactions is modest. However, chemical en- 
ergy from the latent heat of fusion and gravitational en- 
ergy from sedimentation can be important. Understand- 
ing these additional energy sources may allow more ac- 
curate stellar dates, see for example [2J, and provide ad- 
ditional information on the interior composition of WDs 

The interior of a WD is a coulomb plasma of ions and 
a degenerate electron gas. As the star cools this plasma 
crystallizes. Winget et al. recently observed effects from 
the latent heat of crystallization on the luminosity func- 
tion of WDs in the globular cluster NGC 6397 0]. Winget 
et al. 's observations constrain the melting temperature of 
the carbon and oxygen mixtures expected in these WD 
cores. This temperature depends on the ratio of carbon 
to oxygen. Therefore observations of crystallization may 
provide information on WD composition. 

The ratio of carbon to oxygen in WD stars is very 
interesting. It depends on the rate for the reaction 
12 C(a,7) 16 0. Despite a great deal of effort, see for ex- 
ample |5j , the stellar rate for this reaction remains one of 
the most important unsettled rates left in Nuclear Astro- 
physics [BJ. Furthermore, the ratio of carbon to oxygen in 
massive stars is important for their subsequent evolution 
and nucleosynthesis [7]. Therefore, a measurement of the 
carbon to oxygen ratio in a WD could be very important. 

To determine the C/O ratio from observations of the 
melting temperature, one needs the phase diagram for 
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carbon and oxygen mixtures. Recently, we determined 
this phase diagram using molecular dynamics simulations 
and concluded that, if the melting temperature is close 
to that for pure carbon, then the central oxygen fraction 
by mass of the carbon/oxygen WDs in NGC 6397 is less 
than about 64 % [3] . Note that the interior of WDs can 
also be probed with astro-seismology, see for example ref. 

Sedimentation of 22 Ne could provide an additional en- 
ergy source during WD cooling [TU] . Much of the car- 
bon, nitrogen, and oxygen, originally present in the star, 
is converted by nuclear reactions into 22 Ne. This neu- 
tron rich isotope, with 12 neutrons and 10 protons, has a 
larger mass to charge ratio than 12 C or 16 0. As a result 
it will sink in the strong gravitational field of the star 
and release gravitational energy. Recently, Garcia-Berro 
et al. [IT] studied the effects of 22 Ne on WD evolution. 
Sedimentation could be most important in very metal 
rich stars that have more 22 Ne. Both the release of latent 
heat from crystallization and gravitational energy from 
sedimentation can delay WD cooling. Garcia-Berro et al. 
[2] can explain the long 8 Gyr age for the metal rich open 
star cluster NGC 6791 by including both crystallization 
and sedimentation. Furthermore, it may be possible to 
separate these two effects by comparing observations of 
metal rich systems such as NGC 6791 to less metal rich 
systems such as NGC 6397 where sedimentation should 
be unimportant. 

Sedimentation depends on the diffusion constant D for 
22 Ne ions in strongly coupled plasma mixtures. There 
have been previous calculations of D, starting with the 
MD simulations of Hansen et al. for the one compo- 
nent plasma (OCP) [12]. The one component plasma 
consists of ions, with pure coulomb interactions, and an 
inert neutralizing background charge density. Diffusion 
in the OCP in a strong magnetic field was considered by 
Bernu [13]. Hansen et al. have also calculated diffusion 
for binary mixtures [14] . 
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Diffusion for a Yukawa fluid has been simulated by 
Robbins et al. [15] and Ohta et al. [16] , In a Yukawa 
fluid ions interact via a screened coulomb potential (r) , 



Vij(r) 



Z,Z ie 2 



-r/X 
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for two ions with charges Zi and Zj, that are separated 
by a distance r. The OCP is equivalent to a Yukawa 
fluid, where all of the ions have the same charge and the 
screening length A is very large. 

The motion of carbon, oxygen, and neon ions in a WD 
is largely classical because of their large masses. How- 
ever, at great densities there could be some quantum 
corrections that might increase D. These have been es- 
timated by Daligault and Murillo [T7], in a semiclassical 
model, and found to be very small. Instead, these authors 
argue that accurate classical simulations of diffusion co- 
efficients in strongly coupled mixtures with impurities 
would be useful for WD cooling. These diffusion coeffi- 
cients would also be useful to describe sedimentation in 
neutron stars, see for example [IS] . 

In this paper, we present classical MD simulations of 
carbon, oxygen, and neon plasma mixtures in order to 
determine diffusion coefficients D. Our results are more 
accurate than previous work because we explicitly simu- 
late mixtures with realistic WD compositions. In Section 
Il]we describe our MD formalism and present results for 
diffusion coefficents in Section llTll We conclude in Section 

w 



II. FORMALISM 

We describe our MD simulation formalism. This is 
similar to what we used earlier to simulate the carbon/ 
oxygen phase diagram in ref. [3]. We consider a three 
component mixture of carbon ( 12 C), oxygen ( 16 0), and 
neon ( 22 Ne), were the neon is assumed to have a small 
concentration. A star with near solar metallicity, that 
has most of its original carbon, nitrogen, and oxygen 
converted into 22 Ne, might have of order 2% 22 Ne. The 
ratio of carbon to oxygen in the WD core depends on the 
rates for the 4 He(2a,7) 12 C and 12 C(a,7) 16 reactions 
and is expected to be near one to one, see for example 
[31 [20] . Therefore in this paper we present simulations 
using a mixture of 49% 12 C, 49% 16 0, and 2% 22 Ne, by 
number. We expect our results for the Ne self-diffusion 
coefficient to be nearly independent of 22 Ne concentra- 
tion, as long as it is small. Likewise we do not expect 
much dependence of Dno on the ratio of carbon to oxy- 
gen. In addition to this three component mixture, we 
also present results for a one component system of pure 
oxygen for comparison. 

The ions are assumed to interact via screened Yukawa 
interactions, see Eq. [TJ The Thomas Fermi screen- 
ing length A, for cold relativistic electrons, is A -1 = 
2a 1 / 2 kp/'K 1 / 2 where the electron Fermi momentum kp 
is kp = (37r 2 n e ) 1 / 3 and a is the fine structure constant. 



The electron density n e is equal to the ion charge density, 
n e = (Z)n, where n is the ion density and (Z) is the av- 
erage charge. Our simulations are classical and we have 
neglected the electron mass (extreme relativistic limit). 
This is to be consistent with our previous work on neu- 
tron stars. However, the electron mass is important at 
the lower densities in WD and this may change our re- 
sults slightly [21] . Also quantum effects could play some 
role at high densities [221, E3- For relativistic electrons, 
the ratio of A to the ion sphere radius a, 
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depends only on the average charge (Z). For our three 
component mixture we use A = 2.816a, while for pure 
oxygen, we use A = 2.703a. For nonrelativistic electrons 
A/a can be somewhat smaller. In Sec. 



Ill we find that 



our results for diffusion constants D are insensitive to a 
fifty percent decrease in A. 

The simulations can be characterized by an average 
coulomb parameter T, 



(Z 5 / 3 )e 2 
a e T 



(3) 



Here (Z 5 / 3 ) is an average over the ion charges, T is the 
temperature, and the electron sphere radius a e is a e — 
(3/47rn e ) 1 / 3 with n e = (Z)n the electron density. The 
pure system freezes near T = 175 [21] while the C/O/Ne 
mixture is expected to freeze at a somewhat higher T 
[31 [21] . Note that these values of T may depend slightly 
on A [23(25]. 

Time can be measured in our simulations in units of 
one over the plasma frequency uj v . Long wavelength fluc- 
tuations in the charge density can undergo oscillations at 
the plasma frequency. This depends on the ion charge Z 
and mass M. For mixtures we define a hydrodynamical 
plasma frequency Co p from the simple averages of Z and 
M, 
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Note that other choices for the average over composition 
in Eq. [4] are possible. However, they are expected to give 
very similar results for the average plasma frequency. 

Self diffusion constants Z)j, for ions of type i, are cal- 
culated from the velocity autocorrelation function Z l (t), 



Z\t) 



(vi(to)-vj(to)) 



(5) 



where the average is over all ions j of a given type and 
over initial times to- The velocity of the jth ion at time 
t is Vj(t). The diffusion constant is calculated from the 
time integral of Z l {t), 



Di = 



T 

M. 



dtZ\t). 



(6) 
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We find that Z % (t) is small by a maximum time of t max = 
220/a)p. 

We start our MD simulations from a random config- 
uration at a relatively high temperature so that Y ~ 1. 
We evolve the system in time using the simple velocity 
Verlet algorithm [26] with a time step At = l/18w p . We 
use periodic boundary conditions. We do not use a cut- 
off in the interaction at large distances and evaluate the 
force on a given ion by summing over all of the other ions 
in the simulation. We test for finite size effects by per- 
forming simulations for N = 3456, 8192, and 27648 ions. 
Our largest simulation volume is large enough so that 
one half of the box length L is much larger than the elec- 
tron screening length A, L/2 = 8. 65 A. We approximately 
maintain the system at a given temperature by simply 
rescaling the velocities every 10 time steps. Typically we 
equilibrate the system by evolving for a time 2200/w p . 
Then we take data for an additional time of 2200/w p 
by writing the velocities to disk every l/9cD p (two time 
steps). Next we rescale the velocities to a lower temper- 
ature and repeat the equilibration and data taking steps 
for the next higher val ue o f Y. We present results for 
Z l {t) and Di in Section [Till 
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FIG. 1: (Color on line) Velocity autocorrelation function Z l (t) 
versus time t for 22 Ne ions in a carbon, oxygen, and neon 
mixture. The N=8192 ion system is in a fluid state and the 
curves are for values of the coulomb parameter F =9.8 (solid), 
39 (dotted), 98 (dashed), and 244 (dot-dashed). 



III. RESULTS 

We now present results for the velocity autocorrela- 
tion function Z l (t) and diffusion constant Di, first for a 
mixture of carbon, oxygen, and neon and then for a pure 
oxygen system. Fig. [T] shows the velocity autocorrelation 
function for 22 Ne from a simulation with N — 8192 ions. 
The velocity autocorrelation function oscillates with fre- 
quency near uj p and the amplitude of these oscillations 
increases with Y. Note that at Y = 244 the system is still 
in a (possibly metastable) fluid state. 

In Fig. m we show Z l (t) for i = 12 C, 16 0, and 22 Ne 
ions. The difference in Z l for different ions is subtle. We 
see that Z l (t) at t near 5/w p is more negative for 12 C 
than for 22 Ne. Perhaps the lighter carbon ions bounce 
backwards from the confining cage of other ions with a 
more negative velocity than do neon ions. These subtle 
differences in Fig. [2j along with the explicit factor of 
l/Mi in Eq. [6j lead to differences in diffusion constants 
for different ions. 

We choose to scale our diffusion results with Hansen et 
al's simple fit to their original MD results for the diffusion 
constant Dq of a one component plasma |12j . 
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Note that we have generalized ref. [12] results to a multi- 
component system by using the average plasma frequency 
ujp in Eq. [7j In Table |T] and Fig. [3] we present our MD 
simulation results for Di/D . We see that Eq. [7] is in- 
deed a reasonable first approximation to Di and our re- 
sults only differ from this by about 50% or less. Finite 
size effects appear to be modest with 3456 ion simulation 



results being below 8192 or 27648 ion results by only a 
few %. In general there is good agreement between 8192 
and 27648 ion results. However, our results for 22 Ne have 
larger statistical errors than for 12 C or 16 because there 
are fewer 22 Ne ions in our simulations, since its abun- 
dance is assumed to be only 2%. We test the sensitivity 
of our results to the screening length A by running a 
27648 ion C, O, Ne system with a fifty percent smaller 
A. We find at Y — 9.75 that diffusion constants D t only 
increase by two percent or less. 

The dependence of Di on ion species is very interest- 
ing. The 12 C, 16 0, and 22 Ne results in Fig. [3] are nearly 
parallel, for Y > few, and depend on the ion charge Zi 

approximately as Z i . In general, Di depends on both 

the mass and the charge of the ion. However, since the 

ions in our simulations have similar charge to mass ratios 

we do not separate out these two effects and the observed 
2/3 

Z i dependance includes both effects. We can approx- 
imately fit all of the results in Fig. [3] with Y > 5 with a 
simple expression, 



A 



0.53 



(z) 



'(l + 0.22r)exp(-0.135r - 62 ). (8) 



This Eq. is shown as dotted lines in Fig. [3] 

For weak coupling (small Y) the dependance of D j on 
Zi is stronger. In the limit of very weak coupling, Di may 
be related to a mean free path which depends on one over 
an interaction cross section. In Born approximation this 

cross section scales as Zf and gives a Di that depends 

2/3 

more strongly on Zi than Z i . We illustrate this by 
plotting the same data of Fig. [3] but as a function of 1/r 
instead of Y in order to show the low Y data more clearly. 
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TABLE I: Diffusion constants A/A>, for i= 12 C, ie O, and 
22 Ne from simulations with N = 3456, 8192 and 27648 ions, 
for different coulomb parameter F values. The diffusion con- 
stants are scaled with Do = 3ui p a 2 /r 4 ' 3 which is a simple fit 
to the original one component plasma results, see Eq. [7] in 
the text. 

T N 12 C 16 Q 22 Ne 



0.974 



1.392 



1.948 



2.783 



4.871 



9.742 



19.48 



38.97 



64.95 



97.42 



147.6 



3456 
8192 

27648 1.092 0.717 0.460 
3456 
8192 

27648 0.852 0.582 0.449 
3456 0.827 0.580 0.477 
8192 0.827 0.587 0.460 

27648 0.832 0.586 0.453 
3456 
8192 

27648 0.833 0.611 0.500 
3456 0.894 0.697 0.572 
8192 0.881 0.703 0.625 
27648 0.907 0.706 0.609 
3456 1.081 0.888 0.749 
1.058 0.895 0.758 
1.087 0.904 0.788 
1.311 1.121 0.995 
1.343 1.126 1.024 
1.327 1.137 0.984 
1.455 1.251 1.114 
1.464 1.232 1.098 
1.488 1.278 1.119 
1.446 1.216 1.076 
1.420 1.229 1.120 
1.446 1.228 1.100 

1.294 1.077 0.940 
1.297 1.082 0.918 

1.295 1.093 0.979 



8192 
27648 
3456 
8192 
27648 
3456 
8192 
27648 
3456 
8192 
27648 
3456 
8192 
27648 
3456 
8192 

27648 1.034 0.859 0.718 
194.84 3456 

8192 0.779 0.628 0.547 
27648 0.783 0.638 0.533 
221.40 3456 
8192 

27648 0.646 0.516 0.443 
243.55 3456 

8192 0.526 0.408 0.337 
27648 0.544 0.433 0.352 
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FIG. 2: (Color on line) Velocity autocorrelation function Z l (t) 
versus time t for i = 22 Ne (solid cuve), 16 (dotted), and 12 C 
(dashed) ions. The simulation has N = 8192 ions at coulomb 
parameter T = 195.2 and is in a liquid state. 
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FIG. 3: (Color on line) Diffusion constants Di over a simple 
fit to the original MD results for the one component plasma 
Do, see Eq. Ijl versus coulomb parameter F. Results for 12 C 
(solid black line), le O (dashed red), and 22 Ne (dot-dashed 
blue) are presented for a simulation with 27648 ions. Results 
for 8192 ions (squares) are also shown. Finally the dotted 
lines show the simple global fit of Eq. [8] 



We see that Di/D increases for 12 C for small T in a way 
that is not described by Eq. [5J Indeed Fig. [4] shows that 
Di depends more strongly on Zi than Eq. (8ffor small T. 

Finally we also perform MD simulations for a one 
component system of pure 16 in order to compare to 
our multicomponent results. Results for MD simula- 
tions with 8192 and 27648 ions are shown in Table [Til 
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FIG. 4: (Color on line) Diffusion constants Di over Do, see Eq. 
[7| versus the inverse of the coulomb parameter 1/F. Results 
for 12 C (solid black line), 16 (dashed red), and 22 Ne (dot- 
dashed blue) are presented for a simulation with 27648 ions. 
Results for 8192 ions (squares) are also shown. Finally the 
dotted lines show the simple global fit of Eq. [8] This disagrees 
with 12 C and 16 O data for small V 



FIG. 5: (Color on line) Diffusion constant Di over Do, see Eq. 
[7J versus the coulomb parameter F for pure oxygen. Results 
for a simulation with 27648 ions are shown as circles connected 
by a solid black line while the red squares show 8192 ion 
results. Finally the dashed line shows the simple global fit of 
Eq. [8] with {Z)/Zi = 1. 



TABLE II: Diffusion constant Dib /Dq, for a pure le O sys- 
tem from simulations with iV = 8192 and 27648 ions in a 
(possibly metastable) fluid state. The diffusion constant is 
scaled with Do = 3ui p a 2 /F 4 ^ 3 . This is a simple fit to the 
original one component plasma results, see Eq. [7] in the text. 

D/D D/D 
T TV = 8192 N=27648 



2.470 0.680 

6.175 0.831 

12.35 1.035 

24.70 1.237 

49.40 1.307 

82.33 1.194 

123.50 1.000 
187.12 

247.00 .439 



0.828 
1.045 
1.274 
1.330 
1.223 
1.014 
0.703 
0.441 



and Fig. [5j There is good agreement between these two 
simulations and reasonable agreement with Eq. [8] with 
(Z)/Zi = 1. Previous MD simulations for one component 
Yukawa fluids fit D with [H US UHl 



D 



0.0028 



0.00525 ( 



173 



1.154 



(9) 



The ratio of Eq. [9] to Eq. [7] is very similar in form to Fig. 
[5] see Fig. 5 of Ref. [T7]. However, we find a somewhat 
larger amplitude for the deviations of D from Dq ■ 



IV. CONCLUSIONS 

Sedimentation of the neutron rich isotope 22 Ne may 
be an important source of gravitational energy during 
the cooling of white dwarf stars. This depends on the 
diffusion constant for 22 Ne in strongly coupled plasma 
mixtures. We have calculated self-diffusion constants Di 
from molecular dynamics simulations of carbon, oxygen, 
and neon mixtures. 

We find that, for strong coupling (coulomb parameter 
r > few), Di has a modest dependence on the charge 
Zi of the ion species, Di cx Z i ' . However Di depends 
more strongly on Zi for weak coupling (smaller T). Our 
results for both a carbon, oxygen, neon mixture, and for 
a one component plasma can be fit by, 



A 

UJpCL 1 



1.59 



(- 



0.22r 



p4/3 



) exp(- 



0.135r u ' b2 ) 



'(f) 



2/3 



for T > few, see Eqs 
frequency and a is t 



8|7 



(10) 

Here to p is the average plasma 
;hc ion sphere radius. We conclude 
that the self-diffusion constant D^e for 22 Ne in carbon, 
oxygen plasma mixtures is accurately known so that un- 
certainties in Djve should be unimportant for simulations 
of white dwarf cooling. 

All of our diffusion results have been for the liquid 
phase. We are not aware of any results for Di in solids. 
Often Di is arbitrarily assumed to be zero for a solid. In 
future work we will simulate Di in the solid phase, both 
for single component and multicomponent systems. 
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